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We introduce a new class of cellular automata to model 
reaction-diffusion systems in a quantitatively correct way. 
The construction of the CA from the reaction-diffusion equa- 
tion relies on a moving average procedure to implement dif- 
fusion, and a probabilistic table- lookup for the reactive part. 
The applicability of the new CA is demonstrated using the 
Ginzburg-Landau equation. 
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Cellular automata models have been used in many ap- 
plications to model reactive and diffusive systems [|l],D. 
Most uses of cellular automata (CAs) can be classified 
into one of four approaches: (i) Ising-type models of 
phase transitions; (ii) lattice gas models (the lattice gas 
method was initially developed to model hydrodynamic 
flows and has been extended in many directions 

(iii) systematic investigation of the behavior of CAs by 
investigating all rules of a certain class (e.g., all possi- 
ble rules for one-dimensional automata with two states 
and nearest neighbor interaction) [Q; and (iv) qualita- 
tive discrete modelling (including operational use of CAs 
as an alternative to partial differential equations (PDEs) 
1^). These models are generally based on qualitative 
rather than quantitative information about the system 
to be modelled. A CA is constructed which preserves 
the qualitative features deemed most relevant and it is 
then investigated (The lattice gas methods were also de- 
veloped along these lines 0). Existing CA models for 
reaction-diffusion systems |[^~pT|| fall into the category 

(iv) , i.e., they show qualitatively "correct" behavior and 
are restricted to certain reaction-diffusion (R-D) models 
and certain types of phenomena. This is the main crit- 
icism of experimentalists and researchers working with 
partial differential equation models, who search for quan- 
titative predictions. In this letter we describe a class 
of CAs which is suitable for modelling many reaction- 
diffusion systems in a quantitatively correct way. The 
new CAs are operationally more efficient than the reac- 
tive lattice gas methods, which also achieve quantitative 
correctness. We first describe the construction of the 
new class of CAs; then we present the automaton using 
the Ginzburg-Landau equation as an example. 

The main idea behind this class of CAs is careful dis- 
cretization. Space and time are discretized as in normal 



finite difference methods for solving the PDE's. Finite 
difference methods then proceed to solve the resulting 
coupled system of TV x s ordinary differential equations 
{N points in space, s equations in the PDE system) by 
any of a number of numerical methods, operating on 
floating point numbers. The use of floating point num- 
bers on computers implies a discretization of the contin- 
uous variables. The errors introduced by this discretiza- 
tion and the ensuing roundoff errors are often not con- 
sidered explicitly, but assumed to be small because the 
precision is rather high (8 decimal digits for usual float- 
ing point numbers). In contrast, in the CA approach, all 
variables are explicitly discretized into relatively small 
integers. This discretization allows the use of lookup ta- 
bles to replace the evaluation of the nonlinear rate func- 
tions. It is this table lookup, combined with the fact 
that all calculations are performed using integers instead 
of floating point variables, that accounts for an improve- 
ment in speed of orders of magnitude on a conventional 
multi-purpose computer. The undesirable effects of dis- 
cretization are overcome by using probabilistic rules for 
the updating of the CA. 

The state of the CA is given by a regular array of con- 
centration vectors y residing on a d-dimensional lattice. 
Each y(r) is a s-vector of integers (s is the number of 
reactive species). For reasons of efficiency, and to ful- 
fill the finiteness condition of the definition of cellular 
automata, each component yi(r) can only take integer 
values between and b^, where the b^'s can be different 
for each species i. The position index r is a d-dimensional 
vector in the CA lattice. For cubic lattices, r is a d- vector 
of integers. 

The central operation of the automaton consists of cal- 
culating the sum 

y,(r)= ^ y,(r + r') (1) 

r'eNi 

of the concentrations in some neighborhood Ni. The 
neighborhoods can be different for each species i. A 
neighborhood is specified as a set of displacement vec- 
tors, e.g. in two dimensions 

^5-.,a. = {(°),(;,),(;),(V),(-°r)} (2) 

or 

-^1— square — 
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iV5-sta.u{(;),(_\), (-),(:;)}. (3) 

For diffusive systems, symmetry requires that 

r' eN ^ -r' eN (4) 

and if, as is usually the case in reaction-diffusion systems, 
isotropy is required, then 

iri,---,rdf e (r„,,---,r„J^ £ N (5) 

for all index permutations (ai, • • • , ad). For some neigh- 
borhoods, the summing operation can be executed in an 
extremely efficient way by using moving averages: in one 
dimension, the sum of all 2k + 1 cells centered around cell 
r can be computed from the corresponding sum centered 
around cell r — 1 with just one addition and subtraction: 

Yiir) = ^ yt{r + r') 

\r'\<k 

= yt{r -k)+ y,(r - fc + 1) + • • • y,(r + k) 
= -yt{r -k-l)+ y,{r - k - 1) + y,{r - k) + 

•••+y,(r + fc-l)+y,(r + fc) 
= -y^{r -k-l)+ y,(r - 1) + y,(r + k). (6) 

Using this relationship recursively, and using the fact 
that the sum over a square neighborhood in d dimen- 
sions can be constructed as the convolution of such one- 
dimensional sums applied in each dimension in turn, one 
obtains an algorithm to compute the sum of all cells in a 
square (cubic) neighborhood with only 2d additions per 
cell [H . In a multispecies model several variables are nec- 
essary, but they can be packed into one computer word. 
In this case the averaging operation can be performed 
on several species at once if the diffusion coefficients are 
equal. 

In the following we use the normalized values Xi(r) = 
y^{r)/h^ and Xj(r) = yi{r)/{hi\Ni\), which are always 
between zero and one. The resulting fields Xi(r) are then 
the local averages of the Xi(r). The averaging has the 
effect of diffusion. This can be seen from a Taylor ex- 
pansion of Xi(r + r') around x,;(r): 
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\Ni\ ^ ^ k\ y dr 
Xi{r) + DiV^Xiir) + ■ 



(7) 



The factors Di can be computed as shown in |8| and are 
easily calculated from (j^) for square neighborhoods with 
radius k: Di = k{k + l)/6. 

The second operation in the cellular automaton is the 
implementation of the reactive processes described by a 
rate law. Given the reaction-diffusion equation 
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*(z) -fPV^z 



(8) 



(where is understood to be a spatial operator acting 
on each component of the vector z separately, and T) is 
a diagonal matrix), we discretize the time derivative to 
obtain 



= z* + At*(z* ) + AtPV^z* 



(9) 



Changing the time and space scales by setting t ~ t* / At 
and r — r*/Ar, and using the variable x for the rescaled 
set gives 



x*+i = X* + At*(x* 



At_ 



PV^x*. 



(10) 



as the equation to be treated by the CA. Let us define 

**(x) =x + Ai*(x). (11) 
From eqns. (0) and (|ll|) 

**(i*) = X* + AV^x* + ■ • • + At*(x* + AV^x* + ■ • •) 



= X* + AV^x* + At*(x*) + 0{At^ 



Then 



**(x*) 



(12) 



(13) 



is consistently first order accurate in time and within this 
limit, eq. (|l^) can be validly identified with eq. ( p^ to 
describe the evolution of the system. The identification 
yields 



At 

Ar^ 
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(14) 



which defines the space scale. As y is the result of the 
diffusion step, the average output of the CA reaction- 
diffusion process should therefore be given by 



hAN,, 



(15) 



for species j. |^ 

A simplistic discretization of the rate law may pro- 
duce problems: due to the discretization spurious steady 
states or oscillations can appear. It is at this stage that 
the probabilistic rules come into effect. Given an input 
configuration y*(r), one assigns new values y*+^(r) prob- 
abilistically in such a way that the average result corre- 
sponds to the finite difference approximation to the given 
reaction-diffusion equation, gi. The simplest CA rule for 
the reactive step is to treat each species separately and 
+ 1 with probability (gj mod 1) and 
[gjj otherwise ( [ej is the largest integer smaller 
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^Possibly the function needs to be truncated to conform 
to the condition ** e [0, 1]". 
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than or equal to e) . In this way the average is exactly gj . 
Clearly this method introduces the minimal amount of 
noise for a given set of b^'s. In case higher noise levels are 
desired (e.g., if one wants to evaluate the role of fluctua- 
tions ) , one can choose a different set of output values 
and associated probabilities with the same average gj, 
but different variance. To do so, in general at least three 
different possible outcomes (e.g. [gjj — 1, [g^J , [gjj + 1) 
are necessary. 

To demonstrate the applicability of the new CA we 
show how the Ginzburg-Landau equation can be mapped 
onto the automaton. We consider the PDE p^ , p^ 

^^DV^z + az-b\z\^z, (16) 

which can be viewed (for D real) as a two-species 
reaction-diffusion system by separating real and imagi- 
nary parts as z — X + iy, a — a + ij, and b — P + iS, 

dx 

— = DV'^x + ax-^y+ (-f]x + Sy)(x'^ + y^) (17a) 
ot 

^ = DV^y + ay + jx+ {-(3y - Sx)ix^ + y^) (17b) 

where x, y are space- and time-dependent variables. 
With z = re*"^, the corresponding rate equation can be 
conveniently written as 

^ ar - Pr^ = Pr{a/P -r'^) (18a) 

^=1-Sr'. (18b) 

One can set /3 = 1 by changing the time scale. The stable 
homogeneous solution is z = rgC*^* with = \J aj (i and 
ri = 7 — (5r^ (A steady but unstable solution is z = 
0). Since the equation for r is independent of 0, one 
can transform the solution for = to any given il by 
multiplying z(t) with a factor e"^*, yielding an oscillating 
solution. The full PDE system, eq. ([l7|), also admits 
oscillating and rotating spirals, and other inhomogeneous 
solutions. 

For the CA simulations, the region of interest in con- 
centration space needs to be included in the unit square 
[0, 1]^ by shifting the origin to (0.5, 0.5) and scaling the 
parameters in such a way that < Rg < 0.5 (Here 
Rs = 0.25 and j3 ~ 1). In all CA simulations we use a 
square neighborhood (radius k — 1), i.e., Di — D2 — 1/3, 
and discretization levels bi = b2 = 55 (which gives a 
lookup-table size of 9 Mbytes). 

Study of the homogeneous solutions shows that the 
CA behaves like an explicit finite difference method with 
added noise. The time step must be sufficiently small 
to reproduce the correct solution. The noise, which is 
intrinsic in the automaton and arises from the discretiza- 
tion, has little effect on the amplitude of the solution, 
but it introduces random drifts in the phase. 
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FIG. 1. Spiral wave for At — 2, system size 200^ cells. 
Superimposed on the grayscale plot of x is the fitted Archi- 
median spiral. 

We now turn to more interesting cases by considering 
spiral wave solutions. Without loss of generality we use 
f2 = 0, which eliminates the homogeneous oscillations 
and thereby allows for larger time steps in the computa- 
tion and better visualization. 

We initiate a spiral in a two-dimensional simulation by 
starting with the initial condition 

Rirx,ry, i") = ci^(r, - rS)^ + (r, - r")^ (19a) 

dirxTry, t^) — C2 arctan ^, (19b) 

rx - r'i 

which creates exactly one phase singularity at (f^j'^y)- 
For nonzero S smaller than a critical value, a spiral de- 
velops and rotates steadily after some time. In Figure 
|l| we show the spiral obtained for (5 = 1, At = 2, and 
system size (200Aa;)^. Notice that in spite of the use of 
a square neighborhood for the CA, the spiral is perfectly 
round. While this feature is expected when the CA is 
viewed as a numerical method for solving the PDE, it is 
not universally accepted that CA dynamics can produce 
such an isotropic behavior. 

In order to determine the effect of using large time 
steps, we measure the asymptotic wavelength A by fit- 
ting an Archimedian spiral to the contour a; < 0,?; = 0. 
The measured wavelengths for different time steps At 
are shown in Figu re ^ . The value expected from the the- 
ory described in |14| ] is around A*^ — 81. We observe 
that the expected value is reached for small enough At. 
Even for bigger At the deviations are relatively small 
and approximately linear in At. When the parameter S 
becomes bigger than some critical value (estimated to be 
5c = P 1.397 ... in |Q), the spiral wave solution becomes 
unstable. 
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FIG. 2. Spiral wavelength as a function of the size of the 
time step At (logarithmic scale) . The data can be reasonably 
approximated as A « 80.6 + 1.7Ai (solid curve). The size of 
the system is selected to be (6Ao)^ ~ 490^ in non-dimensional 
units, corresponding to 490/\/3AI cells and the spiral fit is 
performed at t=4500. 

As an example of greater coniplexity, Figure ^ shows 
the simulation of a large system (SOO^cells) initialized in 
the state z w 0: Under such conditions many interacting 
spirals develop. Indeed, the initial state is unstable, and 
because of the intrinsic (low level) noise, different regions 
depart from this unstable state with different phase val- 
ues <j). This situation automatically creates many phase 
singularities (points with r = 0, surrounded by points 
with all values of (/>), which then develop into spirals. 
These phase singularities can merge and move as they 
arc influenced by each other |^ . 

Such large simulations are made possible by the speed 
advantage that the CA offers over other numerical meth- 
ods for solving PDEs. We found that on a NeXT work- 
station the CA runs about 5 times faster than an Euler 
integration using the same time step. This advantage is 
due to the use of integer arithmetic, table lookup, and 
the fact that the diffusion calculation is faster in the CA 
than in finite difference methods. The speed advantage 
is much more pronounced when the nonlinear rate func- 
tions are more complicated, as the table lookup operation 
does not slow down. 

In conclusion, we have constructed a class of cellu- 
lar automata that can be used to model many reaction- 
diffusion systems with quantitatively correct results. The 
solutions are very isotropic despite the discreteness and 
the anisotropic nature of the CA. The adverse effects of 
discretization are overcome by the use of probabilistic 
rules. We have demonstrated the applicability of the au- 
tomaton for the Ginzburg-Landau equation. We have 
also successfully applied the method to more than ten 
different reaction-diffusion models, some of which will be 
reported in a forthcoming paper. 
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FIG. 3. Collection of many interacting spirals in a large 
domain (500^cells, 5 = 1, At = 2, t=20000) 
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